home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / src / dassl.cc < prev    next >
C/C++ Source or Header  |  1996-11-03  |  8KB  |  368 lines
  1. /*
  2.  
  3. Copyright (C) 1996 John W. Eaton
  4.  
  5. This file is part of Octave.
  6.  
  7. Octave is free software; you can redistribute it and/or modify it
  8. under the terms of the GNU General Public License as published by the
  9. Free Software Foundation; either version 2, or (at your option) any
  10. later version.
  11.  
  12. Octave is distributed in the hope that it will be useful, but WITHOUT
  13. ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  14. FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  15. for more details.
  16.  
  17. You should have received a copy of the GNU General Public License
  18. along with Octave; see the file COPYING.  If not, write to the Free
  19. Software Foundation, 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
  20.  
  21. */
  22.  
  23. #ifdef HAVE_CONFIG_H
  24. #include <config.h>
  25. #endif
  26.  
  27. #include <string>
  28.  
  29. #include <iostream.h>
  30.  
  31. #include "DASSL.h"
  32.  
  33. #include "defun-dld.h"
  34. #include "error.h"
  35. #include "gripes.h"
  36. #include "help.h"
  37. #include "oct-obj.h"
  38. #include "pager.h"
  39. #include "pt-fvc.h"
  40. #include "utils.h"
  41. #include "variables.h"
  42.  
  43. // Global pointer for user defined function required by dassl.
  44. static tree_fvc *dassl_fcn;
  45.  
  46. static DASSL_options dassl_opts;
  47.  
  48. ColumnVector
  49. dassl_user_function (const ColumnVector& x, const ColumnVector& xdot, double t)
  50. {
  51.   ColumnVector retval;
  52.  
  53.   int nstates = x.capacity ();
  54.  
  55.   assert (nstates == xdot.capacity ());
  56.  
  57.   octave_value_list args;
  58.   args(2) = t;
  59.  
  60.   if (nstates > 1)
  61.     {
  62.       Matrix m1 (nstates, 1);
  63.       Matrix m2 (nstates, 1);
  64.       for (int i = 0; i < nstates; i++)
  65.     {
  66.       m1 (i, 0) = x (i);
  67.       m2 (i, 0) = xdot (i);
  68.     }
  69.       octave_value state (m1);
  70.       octave_value deriv (m2);
  71.       args(1) = deriv;
  72.       args(0) = state;
  73.     }
  74.   else
  75.     {
  76.       double d1 = x (0);
  77.       double d2 = xdot (0);
  78.       octave_value state (d1);
  79.       octave_value deriv (d2);
  80.       args(1) = deriv;
  81.       args(0) = state;
  82.     }
  83.  
  84.   if (dassl_fcn)
  85.     {
  86.       octave_value_list tmp = dassl_fcn->eval (0, 1, args);
  87.  
  88.       if (error_state)
  89.     {
  90.       gripe_user_supplied_eval ("dassl");
  91.       return retval;
  92.     }
  93.  
  94.       if (tmp.length () > 0 && tmp(0).is_defined ())
  95.     {
  96.       retval = tmp(0).vector_value ();
  97.  
  98.       if (error_state || retval.length () == 0)
  99.         gripe_user_supplied_eval ("dassl");
  100.     }
  101.       else
  102.     gripe_user_supplied_eval ("dassl");
  103.     }
  104.  
  105.   return retval;
  106. }
  107.  
  108. DEFUN_DLD (dassl, args, ,
  109.   "dassl (\"function_name\", x_0, xdot_0, t_out)\n\
  110. dassl (F, X_0, XDOT_0, T_OUT, T_CRIT)\n\
  111. \n\
  112. The first argument is the name of the function to call to\n\
  113. compute the vector of residuals.  It must have the form\n\
  114. \n\
  115.   res = f (x, xdot, t)\n\
  116. \n\
  117. where x, xdot, and res are vectors, and t is a scalar.")
  118. {
  119.   octave_value_list retval;
  120.  
  121.   int nargin = args.length ();
  122.  
  123.   if (nargin < 4 || nargin > 5)
  124.     {
  125.       print_usage ("dassl");
  126.       return retval;
  127.     }
  128.  
  129.   dassl_fcn = is_valid_function (args(0), "dassl", 1);
  130.   if (! dassl_fcn)
  131.     return retval;
  132.  
  133.   ColumnVector state = args(1).vector_value ();
  134.  
  135.   if (error_state)
  136.     {
  137.       error ("dassl: expecting state vector as second argument");
  138.       return retval;
  139.     }
  140.  
  141.   ColumnVector deriv = args(2).vector_value ();
  142.  
  143.   if (error_state)
  144.     {
  145.       error ("dassl: expecting derivative vector as third argument");
  146.       return retval;
  147.     }
  148.  
  149.   ColumnVector out_times = args(3).vector_value ();
  150.  
  151.   if (error_state)
  152.     {
  153.       error ("dassl: expecting output time vector as fourth argument");
  154.       return retval;
  155.     }
  156.  
  157.   ColumnVector crit_times;
  158.   int crit_times_set = 0;
  159.   if (nargin > 4)
  160.     {
  161.       crit_times = args(4).vector_value ();
  162.  
  163.       if (error_state)
  164.     {
  165.       error ("dassl: expecting critical time vector as fifth argument");
  166.       return retval;
  167.     }
  168.  
  169.       crit_times_set = 1;
  170.     }
  171.  
  172.   if (state.capacity () != deriv.capacity ())
  173.     {
  174.       error ("dassl: x and xdot must have the same size");
  175.       return retval;
  176.     }
  177.  
  178.   double tzero = out_times (0);
  179.  
  180.   DAEFunc func (dassl_user_function);
  181.   DASSL dae (state, deriv, tzero, func);
  182.   dae.copy (dassl_opts);
  183.  
  184.   Matrix output;
  185.   Matrix deriv_output;
  186.  
  187.   if (crit_times_set)
  188.     output = dae.integrate (out_times, deriv_output, crit_times);
  189.   else
  190.     output = dae.integrate (out_times, deriv_output);
  191.  
  192.   if (! error_state)
  193.     {
  194.       retval.resize (2);
  195.  
  196.       retval(0) = output;
  197.       retval(1) = deriv_output;
  198.     }
  199.  
  200.   return retval;
  201. }
  202.  
  203. typedef void (DASSL_options::*d_set_opt_mf) (double);
  204. typedef double (DASSL_options::*d_get_opt_mf) (void);
  205.  
  206. #define MAX_TOKENS 3
  207.  
  208. struct DASSL_OPTIONS
  209. {
  210.   const char *keyword;
  211.   const char *kw_tok[MAX_TOKENS + 1];
  212.   int min_len[MAX_TOKENS + 1];
  213.   int min_toks_to_match;
  214.   d_set_opt_mf d_set_fcn;
  215.   d_get_opt_mf d_get_fcn;
  216. };
  217.  
  218. static DASSL_OPTIONS dassl_option_table [] =
  219. {
  220.   { "absolute tolerance",
  221.     { "absolute", "tolerance", 0, 0, },
  222.     { 1, 0, 0, 0, }, 1,
  223.     DASSL_options::set_absolute_tolerance,
  224.     DASSL_options::absolute_tolerance, },
  225.  
  226.   { "initial step size",
  227.     { "initial", "step", "size", 0, },
  228.     { 1, 0, 0, 0, }, 1,
  229.     DASSL_options::set_initial_step_size,
  230.     DASSL_options::initial_step_size, },
  231.  
  232.   { "maximum step size",
  233.     { "maximum", "step", "size", 0, },
  234.     { 2, 0, 0, 0, }, 1,
  235.     DASSL_options::set_maximum_step_size,
  236.     DASSL_options::maximum_step_size, },
  237.  
  238.   { "relative tolerance",
  239.     { "relative", "tolerance", 0, 0, },
  240.     { 1, 0, 0, 0, }, 1,
  241.     DASSL_options::set_relative_tolerance,
  242.     DASSL_options::relative_tolerance, },
  243.  
  244.   { 0,
  245.     { 0, 0, 0, 0, },
  246.     { 0, 0, 0, 0, }, 0,
  247.     0, 0, },
  248. };
  249.  
  250. static void
  251. print_dassl_option_list (ostream& os)
  252. {
  253.   print_usage ("dassl_options", 1);
  254.  
  255.   os << "\n"
  256.      << "Options for dassl include:\n\n"
  257.      << "  keyword                                  value\n"
  258.      << "  -------                                  -----\n\n";
  259.  
  260.   DASSL_OPTIONS *list = dassl_option_table;
  261.  
  262.   const char *keyword;
  263.   while ((keyword = list->keyword) != 0)
  264.     {
  265.       os.form ("  %-40s ", keyword);
  266.  
  267.       double val = (dassl_opts.*list->d_get_fcn) ();
  268.       if (val < 0.0)
  269.     os << "computed automatically";
  270.       else
  271.     os << val;
  272.  
  273.       os << "\n";
  274.       list++;
  275.     }
  276.  
  277.   os << "\n";
  278. }
  279.  
  280. static void
  281. set_dassl_option (const string& keyword, double val)
  282. {
  283.   DASSL_OPTIONS *list = dassl_option_table;
  284.  
  285.   while (list->keyword != 0)
  286.     {
  287.       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
  288.                 list->min_toks_to_match, MAX_TOKENS))
  289.     {
  290.       (dassl_opts.*list->d_set_fcn) (val);
  291.  
  292.       return;
  293.     }
  294.       list++;
  295.     }
  296.  
  297.   warning ("dassl_options: no match for `%s'", keyword.c_str ());
  298. }
  299.  
  300. static octave_value_list
  301. show_dassl_option (const string& keyword)
  302. {
  303.   octave_value_list retval;
  304.  
  305.   DASSL_OPTIONS *list = dassl_option_table;
  306.  
  307.   while (list->keyword != 0)
  308.     {
  309.       if (keyword_almost_match (list->kw_tok, list->min_len, keyword,
  310.                 list->min_toks_to_match, MAX_TOKENS))
  311.     {
  312.       return (dassl_opts.*list->d_get_fcn) ();
  313.     }
  314.       list++;
  315.     }
  316.  
  317.   warning ("dassl_options: no match for `%s'", keyword.c_str ());
  318.  
  319.   return retval;
  320. }
  321.  
  322. DEFUN_DLD (dassl_options, args, ,
  323.   "dassl_options (KEYWORD, VALUE)\n\
  324. \n\
  325. Set or show options for dassl.  Keywords may be abbreviated\n\
  326. to the shortest match.")
  327. {
  328.   octave_value_list retval;
  329.  
  330.   int nargin = args.length ();
  331.  
  332.   if (nargin == 0)
  333.     {
  334.       print_dassl_option_list (octave_stdout);
  335.       return retval;
  336.     }
  337.   else if (nargin == 1 || nargin == 2)
  338.     {
  339.       string keyword = args(0).string_value ();
  340.  
  341.       if (! error_state)
  342.     {
  343.       if (nargin == 1)
  344.         return show_dassl_option (keyword);
  345.       else
  346.         {
  347.           double val = args(1).double_value ();
  348.  
  349.           if (! error_state)
  350.         {
  351.           set_dassl_option (keyword, val);
  352.           return retval;
  353.         }
  354.         }
  355.     }
  356.     }
  357.  
  358.   print_usage ("dassl_options");
  359.  
  360.   return retval;
  361. }
  362.  
  363. /*
  364. ;;; Local Variables: ***
  365. ;;; mode: C++ ***
  366. ;;; End: ***
  367. */
  368.